|NASA  Contractor  Report  194961 
ICASE  Report  No.  94-68 


¥ 


ICASE 


A  STABLE  PENALTY  METHOD  FOR  THE 
COMPRESSIBLE  NAVIER-STOKES  EQUATIONS. 
L  OPEN  BOUNDARY  CONDITIONS 


J.  S.  Hesthaven 
D.  Gottlieb 


Contract  NASl-19480 
July  1994 


94-31018 


Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 

Hampton,  VA  2368 1  -0001  SnC  QDALCT7  HTSPECrao  d 


Operated  by  Universities  Space  Research  Association 

94  9  2  8  0  4  0 


A  Stable  Penalty  Method  for  the  Compressible 
Navier-Stokes  Equations. 

I.  Open  Boundary  Conditions. 


J.  S.  Hesthaven  and  D.  Gottlieb  ^ 

Division  of  Applied  Mathematics 
Brown  University 
Providence,  RI  02912,  USA. 

Abstract 

The  purpose  of  this  paper  is  to  present  asymptotically  stable  open  boundary  con¬ 
ditions  for  the  numerical  approximation  of  the  compressible  Navier-Stokes  equations 
in  three  spatial  dimensions.  The  treatment  uses  the  conservation  form  of  the  Navier- 
Stokes  equations  and  utilizes  linearization  and  localization  at  the  boundaries  based  on 
these  variables. 

The  proposed  boundary  conditions  are  applied  through  a  penalty  procedure,  thus 
ensuring  correct  behavior  of  the  scheme  as  the  Reynolds  number  tends  to  infinity.  The 
versality  of  this  method  is  demonstrated  for  the  problem  of  a  compressible  flow  past 
a  circular  cylinder. 
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1.  Introduction.  In  the  present  paper,  we  discuss  boundary  conditions  for  dis¬ 
sipative,  wave  dominated  problems,  exemplified  by  Burgers  equation  and  the  three- 
dimensional,  compressible  Navier-Stokes  equations  given  in  conservation  form.  The 
emphasis  is  on  deriving  open  boundary  conditions  ensuring  the  continuous  problem 
to  be  well-posed  and  on  devising  semi-discrete  schemes  for  imposing  these  conditions, 
which  can  be  proven  asymptotically  stable.  The  boundary  conditions  and  the  semi¬ 
discrete  scheme  are  valid  even  in  the  limit  of  infinite  Reynolds  number. 

When  addressing  exterior,  wave-dominated,  dissipative  problems,  one  is  often 
forced  to  introduce  an  artificial  boundary  for  computational  reasons.  This  introduces 
the  well  known  problem  of  specifying  appropriate  boundary  conditions  at  the  artificial 
open  boundary.  For  purely  hyperbolic  problems,  it  is  well  known  that  enforcing  these 
boundary  conditions  through  the  characteristic  variables  leads  to  a  stable  approxi¬ 
mation.  However,  for  dissipative  wave  problems  the  procedure  is  considerably  more 


complicated. 

Naturally,  we  must  require  the  boundary  conditions  to  lead  to  a  well-posed  con¬ 
tinuous  problem.  For  wave  problems  of  dissipative  type,  the  problem  must,  in  order  to 
be  compatible  with  weak  boundary  layers,  remain  well-posed  even  in  the  limit  where 
the  dissipation  vanishes  and  the  problem  becomes  purely  hyperbolic.  In  addition  to 
this,  we  wish  the  discrete  approximation  of  the  problem  to  be  asymptotically  stable, 
and  that  the  boundary  conditions  are  easily  implemented. 

For  general  non-linear  problems  the  issues  of  well-posedness  zmd  asymptotic  sta¬ 
bility  are  very  complicated,  and  for  most  problems  relatively  little  is  known.  However, 
as  discussed  by  Kreiss  and  Lorenz  [1],  we  may,  for  a  large  class  of  operators,  simplify 
the  problem  significantly  if  the  solutions  are  smooth.  It  was  shown  that  in  this  case 
it  is  sufficient  to  consider  the  questions  of  well-posedness  and  asymptotic  stability  for 
the  linearized,  constant  coefficient  version  of  the  full  problem. 

The  energy  method  is  applied  to  the  linearized,  constant  coefficient  version  of  the 
continuous  problem  in  order  to  obtain  energy  inequalities  which  bound  the  temporal 
growth  of  the  solutions  to  the  initial-boundary  value  problem.  This  technique  allows 


for  handling  such  complex  problems  as  the  Navier-Stokes  equations  and  is  in  general 
applicable  to  symmetrizable  problems  [2]. 

The  usual  way  to  enforce  the  boundary  conditions  in  the  numerical  scheme,  once 
their  proper  form  for  the  continuous  problem  is  known,  is  to  solve  the  equation  in  the 
interior  of  the  computational  domain,  and  then  enforce  the  boundary  conditions  at  the 
boundary  points.  However,  this  approach  does  not  take  into  account  the  fact  that  the 
equation  should  be  obeyed  arbitrarily  close  to  the  open  boundary.  To  circumvent  tWs 
problem,  Funaro  and  Gottlieb  [3,  4],  and  Carpenter  et  al.  [5]  developed  the  penalty 
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method  which  enforces  the  boundary  conditions,  as  well  as  taking  into  account  the 
equation  at  the  boundary.  They  showed  asymptotic  stability  for  the  scheme  applied  to 
scalar  hyperbolic  equations  and  systems  of  hyperbolic  equations.  Don  and  Gottlieb  [6] 
recently  showed  how  this  idea  can  help  in  applying  the  Legendre  collocation  method 
on  Chebyshev  grids. 

The  proofs  presented  in  this  paper  are  all  done  for  semi-discrete  schemes.  The 
relation  between  the  stability  of  the  semi-discrete  and  the  fuUy  discrete  scheme  was 
recently  discussed  by  Kreiss  and  Wu[7]. 

The  issue  of  well  posed  boundary  conditions  for  the  compressible  Navier-Stokes 
equations  was  previously  considered  by  Gustafsson  and  Sundstrom  [8],  Oliger  and 
Sundstrom  [9],  and  Nordstrom  [10].  They  all  used  the  energy  method  to  derive  bound¬ 
ary  conditions  for  the  linearized,  constant  coefficient  Navier-Stokes  equations  in  the 
primitive  variable  formulation.  Dutt  [11]  introduced  an  entropy  function,  which  al¬ 
lowed  him  to  derive  boundary  conditions  for  the  non-linear  problem,  ensuring  that 
the  solution  remains  bounded  in  an  entropy  norm. 

The  remaining  part  of  this  paper  is  organized  as  follows.  In  Sec.  2  we  review 
some  well  known  results  on  Legendre  polynomials  and  collocation  methods.  Section 
3  discusses  Burgers  equation,  and  boundary  conditions  ensuring  weU-posedness  of 
the  problem  are  derived.  We  continue  by  proposing  an  asymptotically  stable  penalty 
method  through  which  the  boundary  conditions  are  enforced.  This  scheme  ensures  the 
correct  behavior  even  in  the  limit  where  the  problem  becomes  hyperbolic,  and  may  in 
general  be  applied  to  any  non-linear  scalar  equation.  The  penalty  method  for  scalar 
hyperbolic,  parabolic,  and  linear  advection-diffusion  equations  is  briefly  discussed, 
and  the  proposed  scheme  is  evaluated  by  numerical  tests.  The  importance  of  properly 
choosing  the  penalty  parameter  is  addressed  in  Sec.  4,  where  we  discuss  the  effect  of 
the  penalty  method  on  the  CFL  condition  when  using  explicit  Runge-Kutta  methods 
for  time-stepping  linear  problems.  The  results  from  the  linear  anedysis  are  s!  h!  own 
to  carry  over  to  the  non-line 

2.  Legendre  Polynomials  and  Collocation  Methods.  The  schemes  which 
we  ansdyze  in  the  present  paper  are  all  based  on  Legendre  collocation  methods.  This 
choice  is  merely  dictated  by  a  wish  to  obtain  analytical  results,  and  the  methods  ex¬ 
tend  trivially  to  other  collocation  methods  and  even  to  finite  difference/finite  element 
methods. 

The  Legendre  polynomial  of  order  N  is  defined  as 

where  |x|  <  1.  We  will  in  the  following  only  consider  collocation  methods,  where  the 


2 


collocation  points  are  given  as  the  Legendre-Gauss-Lobatto  points,  being  defined  as 
the  roots  of  the  polynomial  (1  -  x^)Pjtf(x).  There  is  no  known  explicit  formula  for 
these  roots. 

Associated  with  the  Gauss- Lobatto  points  is  the  quadrature  formula,  stating  that 
if  f(x)  is  a  polynomial  of  degree  2N  —  1,  then 

(1)  = 

*:=0 

where  Xk  are  the  Legendre-Gauss-Lobatto  collocation  points,  and  the  Gauss-Lobatto 
weights,  cjk,  are  given  as 


/,/«)<« . 


For  further  details  on  the  properties  of  the  Legendre  polynomials,  we  refer  to  [12]. 

In  a  Legendre  collocation  method,  the  function,  f(x),  is  approximated  by  a  grid 
function,  /t  =  /(x^),  where  the  grid  points  are  the  Gauss-Lobatto  collocation  points. 
Thus,  we  construct  a  global  Legendre  interpolant,  /;v,  to  obtain  the  approximation  to 
the  function; 


N 


(lNf)(x)  =  '^fkhk(x)  , 


»=0 


where  the  interpolating  Legendre- Lagrange  polynomials  are  given  as 

^  N(N  +  l)(x-Xk)PN(xk)  ■ 

We  note  that  by  construction 


(•fiv/)(®fc)  =  fk  ■ 

To  seek  equations  for  an  approximate  solution,  (/yv/)(®),  to  a  partial  differential 
equation,  we  need  to  obtain  values  for  the  spatial  derivatives  at  the  collocation  points. 
This  is  done  by  approximating  the  differential  operator  by  a  matrix  operator,  with 
the  matrix  entries  given  as 

Pki  =  h\{xk)  ■ 

For  the  explicit  expression  of  the  entries,  we  refer  to  [13,  14]. 
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aU{-l,t)-0e^ 


dx  L=_i 


=  0  . 


^Uil,t)+6e^ 


dx  L=, 


=  0  , 


When  addressing  the  issue  it  is,  as  discussed  in  the  introduction,  sufficient  to  consider 
the  linearized,  constant  coefficient  version  of  Burgers  equation 


Here  X  =  Uo'is  the  uniform  solution  around  which  we  have  linearized.  Equation  (6)  is 
also  known  as  the  linear  advection-diifusion  equation. 

The  four  real  constants,  a,  /?,  7,  and  S,  in  the  boundary  conditions,  Eq.(4)- 
(5),  may  not  be  chosen  arbitrarily,  since  the  resulting  problem  should  be  well-posed. 
Bounds  yielding  a  sufficient  condition  for  well-posedness  are  given  in  the  following 
Lemma. 


Lemma  3.1.  Equation  (6),  with  boundary  conditions  given  by  Eq.(4)-(5)  is  well- 
posed  if  one  of  the  following  conditions  holds 

(i) :  /3  =  0  ,  =  0. 

(ii) :  P  ^0  ,6  =  0  and  (e  -  A)  -|-  2a//3  >  0. 

(Hi):  P  =  0  ,  ^  ^  0  and  (e  -|-  A)  -|-  2j/6  >  0. 

(iv):  p^Q  ,  6^0  and  2{e-  X)'y/6  +  2(£  -h  A)a//3  -|-  A{a'r)/{P6)  >  A^. 

Proof  Construct  the  energy  integral  as 

Here  we  have  introduced 

{U,V)  =  f  uVdx  ,  {U,U)=\\Uf  . 
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Following  the  similar  analysis  done  in  [15],  we  use  the  following  estimate 


Applying  this,  the  condition  for  well-posedness  becomes 
1  d 


Condition  (i)  implies  that  17(-I)  =  £f(l)  =  0  such  that 


For  condition  (ii)  we  obtain  ff(l)  =  0  and  thus 

yielding  the  condition 


A  +  2->0 


Likewise,  for  condition  (iii)  we  obtain 

showing  that  this  choice  yields  well-posedness.  For  condition  (iv)  we  obtain  the  fol¬ 
lowing  condition 

+  2^)  If'C-l)  +etf(-l)(/(l)  -  i  [e  +  A  +  2|)  U\l)  <  0  . 

This  is  obeyed  if 

e^-(£_A  +  25)(5  +  A  +  2j)<0  , 

implying 


2{e  -  X)'y/6 +  2{e  +  X)a/P  ■^A{a'y)/{P6)>  X"^  .  0 

3.1.  The  Semi-Discrete  Scheme.  Equation  (3)  will  be  solved  using  a  Legendre 
collocation  method  where  the  collocation  points  are  the  Legendre-Gauss-Lobatto  points 
This  involves  finding  an  iVth  degree  polynomial,  u(x,t),  satisfying 


(7) 


du  du  _  d^u 
dt  ^  ^dx  ^  dx'^ 


at  X  =  Xfc  ,  A:  e  [1 . .  .IV  -  1]  , 
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in  the  interior.  The  boundary  points  are  given  by  boundary  conditions  of  Robin  type 

au{xoJ)-fi£—  =gi{t)  , 
du 

'yu{xN,t)+  6e-^  =  92{i)  , 

where  gi{t)  and  ^2(^)  aJ'e  the  boundary  conditions.  The  traditional  method  of  imposing 
the  boundary  conditions  is  to  solve  Eq.(7)  in  the  interior  and  enforce  the  boundary 
conditions  at  the  boundary  points  only.  However,  this  approach  does  not  take  into 
account  the  fact  that  the  equation  must  be  obeyed  arbitrarily  close  to  the  boundary. 
In  addition  to  this,  it  has  proven  difficult  to  implement  Robin  boundary  conditions 
consistently  for  non-linear  problems.  To  overcome  these  problems,  we  follow  the  line 
of  thought  initiated  by  Funaro  and  Gottlieb  [3,  4]  and  propose  a  penalty  method 
for  Burgers  equation  at  the  Legendre-Gauss-Lobatto  collocation  points,  x  =  Xk  ,k  S 
[0, . . . ,  J\r],  as 


du  du 
dt  ^  dx 


where 


-  tiQ~{x)  au{xo,t)-P£—  -gi{t) 


-  T2(?+(x)  -52(0 


(9) 


Q  (®) 


(1  -  ^)Pn{^) 


Q^x) 


(1 -I- i)p;^(x) 


These  two  functions  have  the  property  of  being  zero  at  all  Legendre-Gauss-Lobatto 
collocation  points,  except  at  the  two  endpoints  of  the  domain.  Although  Q~  and 
here  are  defined  as  delta-functions  at  the  boundary,  we  may  equally  well  chose  other 
definitions.  As  shown  by  Don  and  Gottlieb  [6],  this  approach  may  sdso  be  applied  for 
implementing  Legendre  methods  on  Chebyshev  grids. 

We  note  here  that  the  penalty  method  as  given  by  Eq.(8)  combines  the  boundary 
conditions  and  the  governing  equation  into  one  equation.  When  using  the  penalty 
method,  the  boundary  conditions  are  not  exactly  obeyed  at  the  boundary.  However, 
the  method  remains  spectrally  ciccurate,  as  we  will  soon  illustrate.  One  may  also 
observe  that  the  scheme  is  eqtiivalent  to  the  traditional  approach  for  ti  ,  T2  approaching 


infinity. 

In  order  to  obtain  the  energy  inequality,  we  consider  only  homogeneous  boundary 
conditions.  As  discussed  previously  in  [1],  this  is  no  restriction,  since  we  may  always 


6 


introduce  a  variable  transform  such  the  boundary  conditions  become  homogeneous. 
In  the  following  Lemma  we  state  the  bounds  on  tj  and  T2  ensuring  that  the  linearized, 
constant  coefficient  version  of  Eq.(8)  is  asymptotically  stable. 


Lemma  3.2.  Assume  u{x,t)  exists  and  let  ^  and  be  defined  as 

r'i  =  +  2k-  +  eK-  l/2fa;|A|]  , 

+  2k  +  2^k2  +  e/c-  l/2ea;|A|]  , 

where  k  =  wa/b  and 

2 

N(N  +  l)  ’ 

is  the  Legendre  weight  at  the  end-points. 

Then  if 

T-0  <  r,  <  r+^  , 

T's  <  -^2  <  <5  , 

then  the  linearized,  constant  coefficient  version  of  Eq.  (8)  is  asymptotically  stable  and 
the  solution  is  bounded  as 


Proof.  We  start  be  defining  the  discrete,  weighted  scalar  product  as 

N 

(u,v)jv  =  J^u(Xk)v(x/i)Wk  ,(«,u)/v  =  ||u11n  • 
k=0 

and  note  that  since  we  are  using  a  Legendre  collocation  method,  we  have,  through 
Eq.(l),  the  identity 

This  makes  it  straightforward  to  apply  partial  differentiation.  Following  the  results 
stated  previously,  it  is  sufficient  to  obtain  the  energy  estimate  for  the  linearized,  con¬ 
stant  coefficient  version  of  Eq.(8); 

-  £1K1In 

-Tiwtt(-l)[au(-l)  -  fieux{-Vj\  -  T2a>u(l)[7u(l)  +  5euj,(l)]  . 
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Here  the  subscripts  designate  differentiation  and  u  is  the  Legendre  weight  at  the 
endpoints  (£q.(2)).  Using  the  quadrature  rule  allows  for  rewriting 

N-l 

k=\ 

Contrary  to  the  approach  followed  by  Funaro  and  Gottlieb  [3,  4],  we  recast  the  prob¬ 
lem  of  stability  into  an  algebraic  eigenvalue  problem.  For  the  present  problem,  this 
may  seem  an  addition2il  complication.  However,  we  find  that  for  more  complicated 
problems,  this  approach  greatly  simplifies  the  proofs. 

Isolating  the  terms  contributing  to  stability  at  each  boundciry,  we  obtain  two 
conditions  for  asymptotic  stability; 

uIW“u-  <  0  ,  uln+u+  <  0  , 

where  u_  =  [«(-!), Ui(-l)]^,  u+  =  [u(l),ti4l)P  and 

A  -  2aujTi  -£(1  —  /JwTi ) 

—£(1  - /3u;ti)  —2£cj 

-A  -  2'yu;r2  £(1  -  tfwTj) 

£(1-SuJT2)  -2£Uf 

Since  both  matrices  are  symmetric,  the  problem  is  reduced  to  ensuring  that  7f~  and 
Tf**  are  negative,  semi-definite.  The  eigenvalues  of  the  two  matrices  are  found  to  be 

^i,2(W")  =  ^  ±  +  16£(/3^ujhTjf  -  2uj{I3£  ->r  2aa;)T]  -|-  2a;A  -f 

^  ±  +  16£(^W£r|  -  2u}{6£  +  2'yu)r2  -  2a;A  -|-  e)^  , 

where  =  -2X  +  ^£  +  4aujTi  and  =  2A-|-4a;£-|-47u;T2.  It  is  evident  that  negative 
semi- definiteness  is  ensured  if 

P^u^£T^  -  2a;(/?c  -H  2aa/)Ti  -|-  2a)A  -|-  e  <  0 
6^ui^£T2  —  2ui{6£  27a; )t2  —  2a;A  -|-  e  <  0  . 

The  roots  of  the  two  polynomials  are 

T*  =  — ^  -f-  2/c_  i:  2^«i  -H  e/c_  —  l/2ea;A^  , 

rf  =  ^£  +  2k+  ±  2^k\  -I-  £K+  +  l/2ea;A^  , 
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where  k_  =  uaj^  and  We  introduce 


~  +  2k  -  2^ +  SK  -  ll2eu)\\\\  , 

T'a^6  =  ;;^[^  +  2k  +  2^k2  +  ck  -  \l2eu)\\\]  , 

wherp  K  =  (4;a/6.  Since 

-  ^  1^1  11 
-  2aa;  4^a;2  ’ 

for  £  <  1,  this  ensures  C“  >  0  and  >  0. 

Hence,  stability  is  ensured  for 

<  ri.<  , 

T's  <T2<  , 

with  the  solution  satisfying 

\ji\MN  -  S  “x(®fcH  .  D 

S.l.l.  Remarks  on  the  Penalty  Method  for  Linear  Equations.  The  results 
stated  in  Lemma  3.2  may  be  used  to  derive  the  appropriate  penalty  parameter  for 
a  large  class  of  linear  equations.  We  consider  the  general  linear  advection-diffusion 
equation,  Ek[.(6),  with  the  Robin  boundary  conditions  given  in  Eq.(4)'(5).  Solving  this 
problem  by  a  penalty  method,  equivalent  to  that  given  by  Eq.(8),  requires  bounds  on 
the  penalty  parameters  in  order  to  ensure  stability  of  the  scheme. 

In  the  following,  we  wiD  give  these  bounds  for  reference  and  will  return  to  the 
numerical  validation  of  these  results  in  Sec.  4.  Some  of  these  results  may  be  found  in 
[3,  4,  6],  but  aue  here  given  in  a  more  general  framework.  Note  that  w  ~  0{N^). 


(i)  Hyperbolic  Equations  (g  =  0). 

1  A  >  0.  Well-posedness  is  ensured  by  choosing  a  >  0  and  (3  =  "y  ~  6  =  0. 
Thus,  for  this  case  we  will  only  need  bounds  on  . 


^a,0  = 


2ti;a 


<0  =  00 


The  scheme  for  the  hyperbolic  case  is  stable  for 


00  >  Ti  > 


A 

2u;a 
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2.  A  <  0.  Well'posedness  is  ensured  by  choosing  7  >  0  and  a  =  ^  =  0. 

Thus,  for  this  case  we  will  only  need  bounds  on  T2. 


=  00  . 


»  S.o 


The  scheme  for  the  hyperbolic  case  is  stable  for 


00  >  T2  > 


(ii)  Parabolic  Equations  (A  =  0,  £  >  0).  Necessary  and  sufficient  conditions  for 
well-posedness  may  be  obtained  by  choosing  the  four  parameters,  a,  /3,  7,  and  S, 
properly  as  stated  in  Lemma  3.1  [15].  We  only  state  the  results  for  the  bounds  of  tj, 
since  the  results  for  are  equivalent. 

1.  Dirichlet  boundary  condition,  (a  >  0,  /3  =  0). 


Stability  is  ensured  for 


-  _  ^  1  +  _ 

"  4aa;2  ’  '^«.o  “  • 


.  .  £  1 
00  >  Ti  >  - - 2  • 

4a  u>^ 


2.  Neumann  boundary  condition,  (a  =  0,  >  0). 


Stability  is  ensured  for 


3.  Robin  boundary  condition,  (a  >  0,  >  0). 


+  2k  -  +  en]  , 

=  -^[£  +  2k  +  2-v/k2  +  £k]  , 


where  k  =  ojafp.  Stability  is  ensured  for 


>  r,  >  . 

(iii)  Advection-Diffusion  Equations  (A  ^  0,  £  >  0).  Agaun  we  must  ensure  well- 
posedness  by  proper  choice  of  the  four  parameters,  as  given  in  Lemma  3.1.  We  only 
state  the  results  for  the  bounds  of  ri  as  the  results  for  t-z  are  equivalent. 
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1.  DiricUet  boundary  condition  (a  >  0,  /3  =  0). 

|A|  1  el  , 

2aa;'*' 404^2  '  '^«,o  -  oo  . 


Stability  is  ensured  for 


lAl  1  £  1 


2.  Neumann  boundary  condition  (a  =  0,  /3  >  0). 
1 


■  I3u 


f2lA|u;  ,+  __L  , 

^2  >  ^0.0  -  13^  + 


2\X\ut 


y32  • 


Stability  is  ensured  for 


Pcj  ^ 


f2|A|a;  .  ^  .  J_  _ 

yg2  -  ^  -  P(j 


l2\X\u; 


/32  • 

3.  Robin  boundary  conditions,  a  >  0,  /3  >  0.  Results  are  given  in  Lemma  3.2. 


3.2.  Numerical  Tests.  As  we  aim  at  solving  the  full  non-linear  Burgers  equa¬ 
tion,  and  not  the  linearized,  constant  coefficient  version,  we  need  to  validate  the  results 
obtained  from  the  linear  analysis.  We  have  solved  Burgers  equation  using  the  scheme 
given  by  Eq.(8)  and  employing  a  standard  Legendre  collocation  method  [13,  16]. 

Burgers  equation,  Eq.(3),  has  a  rightward  traveling  wave  solution  (see  e.g.  [1])  of 
the  form 


(10)  (7(x,t)  = -otanh  -I- c  ,iG  [-00,00]  ,t>0  , 

where  the  free-stream  values 

lim  ff(i,t)  =  6_oo  ,  lim  U{x,t)  =  boo  , 

X—¥^O0  X^OO 

are  associated  with  the  wave-speed,  c,  and  the  constant,  a  >  0,  as 

b—oo  "1"  boo  b—00  boo 

c  = -  ,  o  = -  . 

2  ’  2 

Following  the  results  in  Lemma  3.1  (condition  (iv):  a  =  A,  /3  =  1,  7  =  0,  5  =  1),  we 
expect  the  non-linear  problem  to  be  weU-posed  for  boundary  conditions  of  the  type 

=M<)  . 

where  A  >  0  is  the  value  around  which  we  have  linearized.  In  the  present  study,  we 
have  used  the  free-stream  value  at  the  inflow,  i.e.  A  =  b-oo- 
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Since  we  know  an  exact  solution,  the  boundary  conditions  may  be  given  exactly 
at  all  times  using  Eq.(lO).  As  initial  condition  we  use 

D"(i,0)  =  -otanh  ^  • 

The  solution  is  time-stepped  using  a  classical  4th-order  Runge-Kutta  method,  where 
the  boundary  conditions  are  imposed  at  the  intermediate  time-levels. 

Using  the  values  of  the  penalty  parameters  given  in  Lemma  3.2  results  in  a  stable 
scheme.  However,  the  CFL-number,  relating  the  maximum  2illowable  time  step  to  the 
spatial  resolution  as 

A.  ^  CFL 

1  2  . 

I  I  mui  '  nun 

will  have  to  be  very  small  in  order  to  ensure  stability.  Here  \U\  signifies  the  maximum 
absolute  value  of  U.  Thus,  with  the  theoretical  value  of  the  penalty  parameter,  the 
proposed  method  compares  unfavorably  with  the  traditional  method,  due  to  severe 
time  step  restrictions.  Fortunately,  the  limits  of  the  penalty  parameters,  in  between 
which  asymptotic  stability  is  ensured,  are  obtained  as  a  result  of  a  conservative  energy 
estimate  and  hence  are  not  very  accurate. 

We  have  used  the  values  of  penalty  parameter  (see  Lemma  3.2)  as; 


These  values  are  found  to  lead  to  a  stable  scheme,  provided  eN^  1.  In  Eq.(3),  e 
plays  the  role  of  an  inverse  Reynolds  number.  The  constraint,  eN^  >  1,  simply  states 
that  increasing  the  Reynolds  number  requires  increased  spatial  resolution,  which  is 
a  natural  restriction.  For  advection  dominated  problems,  stability  is  obtained  by 
increasing  the  penalty  parameters  towards  the  values  stated  in  Lemma  3.2. 

With  these  values  of  the  penalty  parameters,  we  have  been  able  to  perform  the 
simulations  with  a  CFL  number  of  4,  which  is  equivalent  to  what  is  usually  allowed 
when  using  the  traditional  method.  Thus,  by  fine-tuning  the  penalty  parameters  we 
were  able  to  avoid  any  effect  of  the  penalty  method  on  the  CFL-condition.  The  follow¬ 
ing  section  contains  a  study  of  the  effect  of  the  penalty  method  on  the  CF’L-condition 
and  guidelines  for  fine-tuning  the  penalty  parameter  for  practical  applications. 

In  Fig.  1  we  show  the  temporal  evolution  of  the  traveling  wave  solution  when 
using  the  proposed  scheme  as  given  by  Eq.(8).  The  simulation  is  done  with  iV  =  64 
and  e  =  0.1.  We  observe  no  spurious  reflections  from  the  open  boundary  and  the  kink 
is  seen  to  travel  undisturbed  out  of  the  domain.  Table  1  shows  the  error  at  T  =  1.00, 
where  the  kink  has  propagated  half  way  through  the  boundary.  It  is  evident  that 
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Pio.  1 .  Traveling  wave  solution  of  Burgers  equation. 

Tablb  1 

Error  in  the  spectral  simulation  of  Burgers  equation  using  the  penalty  method.  The  maximum 
error  (Lao)  occurs  at  the  boundary. 


N 

I-a 

Xoo 

16 

1.07E-02 

3.26E-02 

32 

7.64E-05 

3.50E-04 

64 

3.36B-09 

2.21E-08 

128 

1.58E-11 

7.62E-11 

the  proposed  scheme  maintains  the  spectral  accuracy.  The  time-step  is  so  small  that 
time-stepping  errors  may  be  neglected. 

4.  CFL-Restrictions  for  the  Penalty  Method.  As  discussed  briefly  in  the 
previous  Section,  choosing  too  large  a  penalty  parameter  results  in  severe  CFL- 
restrictions.  For  this  reason,  it  is  vital  to  understand  how  the  penalty  method  alters  the 
eigenvalue  spectrum  of  the  operators  and  consequently  changes  the  CFL- restriction. 

In  the  present  section  we  will  study  these  effects  for  the  linear  advection  and 
diffusion  operators  for  Legendre  collocation  methods.  For  completeness,  we  will  also 
give  the  results  for  Chebyshev  collocation  methods,  which  are  widely  used  when  solving 
non-linear  problems.  The  analysis  will  consider  both  3rd-  and  4th-order  Runge  Kutta 
methods,  which  are  often  employed  when  addressing  problems  of  the  type  considered 
here.  At  the  end  of  the  section  we  will  compare  the  results  from  our  linear  analysis 
with  simulations  of  the  non-linear  Burgers  equation. 

Consider  now  the  semi-discrete  linear,  constant  coefllcient  problem 
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(11) 


(q)t  =  £wq  Xfc  6  n  ,  «  >  0 
q  =  0  XA;6n,t  =  0  , 

Bsq  =  0  x/b€r,<>0 

where  q  =  (q(xo), . .  .,q(x;v))^i  k  G  [0, £jv  is  the  discrete  appronmation  of 
the  operator  for  the  interior  and  Bn  determines  the  appropriate  discrete  boundary 
conditions.  We  assume  that  the  semi-discrete  approximation  is  a  consistent  approxi¬ 
mation  of  the  continuous  problem.  A  time-differencing  scheme,  where  the  boundary 
conditions  are  enforced  exactly  at  the  boundary  points,  may  then  be  expressed  as 

=  KsiAt,  CnK 
=  0  . 

Here  q"  signifies  the  solution  vector  at  time-step  n.  Thus,  for  strong  stability  we  must 
require 


(AyvCAt,  £jv)|  <  1  . 


However,  employing  the  penalty  method  changes  the  time-stepping  scheme  as 

=KNiAtXN-TBN)q'^ 

and  strong  stability  is  ensured  if 


\KN{At  ,Cn  -tBn)\  <  1  , 


explaining  why  the  Cfi-condition  depends  strongly  on  the  correct  choice  of  the 
penalty-puameter. 

In  the  following  analysis  we  consider  explicit  Runge-Kutta  time  stepping  methods, 
which,  for  time  independent  operators,  may  be  expressed  as 


»=o  *' 

where  p  is  the  order  of  the  scheme.  We  have  for  simplicity  assumed  that  the  boundary 
conditions  are  included  in  the  operators.  Assuming  Cn  =  SnAnS]!^^,  where  |5a^|  and 
I  are  bounded  and  independent  of  N,  strong  stability  of  the  Runge-Kutta  schemes 
is  obtained  if 


Cn)\  =  Sn 


£i(A(AA,)‘  V  =  f;i{A(AN)' 

t=0  t=0 


<  1 
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Tablb  2 


Seeding  constants  for  the  advection  operator.  The  proper  boundary  conditions  are  of  Dirichlet  type 

(D). 


Advection  Operator 

Cl 

C 

o' 

A#0  e  =  0 

4tli  RK 

D  Exact  BC 

■Q 

35 

27 

32 

Penalty  BC 

■9 

17 

10 

11 

Hence,  the  problem  is  reduced  to  finding  the  eigenvalue  spectrum  of  the  operator  Cn 
and  choose  At  accordingly. 

In  the  present  study  we  consider  the  linear  advection- diffusion  operator; 

with  the  Robin  boundary  condition  operators 

The  boundary  conditions  for  the  exact  method  are  enforced  through  the  operator  as 
described  in  [16]. 

In  order  to  compare  time-step  restrictions  as  found  for  the  two  different  ap¬ 
proaches,  we  now  define  the  two  CJFi-like  constants,  Ci  and  Cc,  as 

Af  Ai  ^ 

^  ^  -  Aj\r2  -i-  eN*  ’ 

where  the  subscripts  refer  to  Legendre(£)  and  Chebyshev(C')  operators,  respectively. 
These  constants  are  determined  by  solving  the  eigenvalue  problem  and  calculating  the 
maximum  At  which  ensures  stability  and  supplies  an  upper  bound  on  the  time-step. 

Table  2  and  3  shows  the  calculated  values  of  Cl,  and  Cc  for  the  advection  and  the 
diffusion  operator.  The  results  are  the  same  for  the  full  advection-diffusion  operator 
for  the  diffusion  operator,  provided  >  1,  and  is  therefore  omitted. 

It  is  clear  firom  Table  2  that  using  the  penalty  method  for  enforcing  boundary 
conditions  on  purely  advective  problems  results  in  a  significant  reduction  of  the  maxi¬ 
mum  allowable  time-step.  However,  more  importantly,  Table  3  shows  that  for  problems 
where  the  diffusion  operator  dominates  the  eigenvalue  spectrum,  the  penalty  method 
allows  for  increasing  the  time-step  with  as  much  as  50  %.  The  effect  is  most  pro¬ 
nounced  when  using  a  4th-order  Runge-Kutta  method  for  time-stepping  a  Chebyshev 
collocation  scheme. 

In  order  to  explain  the  results  in  Table  2  and  3,  we  compare  in  Fig.  2  the  spectrum 
of  the  Legendre  collocation  advection  (Fig.  2a)  and  diffusion  (Fig.  2b)  operators  when 
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Table  3 


Scaling  constants  for  the  diffusion  operator.  Results  are  given  for  possible  combinations  of  Dirich- 
let  (D),  Neumann  (N)  and  Robin  (R)  boundary  conditions. 


Diffusion  Opei&toi 

A  =  0  e  >  0 

C 

3id  RK 

L 

4th  RK 

C 

3id  RK 

c 

4th  RK 

D-D/D-N/D-R 

Exact  BC 

99 

109 

53 

58 

Penalty  BC 

81 

123 

56 

84 

N-R 

Exact  BC 

99 

109 

53 

58 

Penalty  BC 

130 

135 

91 

96 

R-R 

Exact  BC 

99 

109 

53 

58 

Penalty  BC 

130 

141 

93 

97 

/ 

•  £actec 

•  PonanyBC 

N>24 

’  > 

- 

FlO.  2.  Eigenvalue  spectrum  ('A  s  Ar  +  ihij  for  the  Legendre  advection  operator  (Ba)  and  the 
Legendre  diffusion  operator  (tb)  as  obtained  by  using  exact  boundary  conditions  (•)  and  the  penalty 
method  (o). 

enforcing  Dirichlet  boundciry  conditions  through  the  exact  method  and  the  penalty 
method. 

For  the  advection  operator  (Fig.  2a)  we  observe  that  the  effect  of  the  penalty 
method  is  to  introduce  an  extreme  complex  conjugate  eigenvalue-pair,  which  domi¬ 
nates  the  spectrum  and  eventually  determines  the  maximum  allowable  time-step.  This 
results  in  the  decreased  Cfi-number  as  observed  in  Table  2. 

The  effect  on  the  diffusion  operator  is  more  complicated  and  depends  strongly 
on  the  value  of  the  penalty  parameter.  As  proved  by  Gottlieb  and  Lustman  [15], 
the  diffusion  operator  with  exact  Robin  boundary  conditions  has  a  real,  negative  and 
distinct  eigenvalue  spectrum.  This  property  is  preserved  if  a  sufficiently  large  value  of 
r  is  used  in  the  penalty  method.  However,  by  decreasing  the  penalty  parameter  the 
two  most  extreme  eigenvalues  split  into  two  pairs  of  complex  conjugate  eigenvalues, 
which  move  towards  the  imaginary  axis,  as  r  is  decreased.  In  Fig.  2b  we  show  the 
eigenvalue  spectrum  for  the  optimal  choice  of  r.  The  important  observation  to  make 
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is  that  moduli  of  these  new  eigenvalues  are  smaller  than  the  original  extreme  negative 
real  eigenvalue.  Additionally,  since  the  dominating  eigenvalue  now  is  complex,  it 
clearly  becomes  advantageous  to  use  the  4th-order  Runge-Kutta  method  due  to  the 
increased  extension  of  the  stability  region  along  !  the  imaginary  axis  as  compared  to 
The  validity  of  this  conclusion  is,  however,  strongly  dependent  on  the  proper 
choice  of  the  penalty  parameter.  The  values  derived  in  the  previous  section  do  indeed 
ensure  asymptotic  stability,  but  with  a  significant  reduction  in  the  maximum  allow¬ 
able  Ci^Z>-number  as  a  result.  Fortunately,  as  mentioned  previously,  the  limits  of  the 
penalty  parameters  are  based  on  a  conservative  energy  estimate  and  are  therefore  not 
very  accurate.  In  the  following  we  give  the  penalty  parameters  used  to  obtain  the 
results  given  in  Table  2  and  3.  These  values  result  in  a  stable  scheme  as  long  as  the 
problem  is  purely  advective  or  eN^  >  1,  and  allows  in  most  cases  for  a  significant 
increase  in  the  time-step. 


(i)  Legendre  Collocation  Methods 

1.  Dirichlet  Boundary  Conditions. 

2.  Neumann  Boundary  Conditions. 

N{N  -I- 1) 
^  ~  8 

3.  Robin  Boundary  Conditions. 


r  = 


(ii)  Chebyshev  Collocation  Methods 

1.  Dirichlet  Boundary  Conditions. 


2.  Neumann  Boundary  Conditions. 


T  = 


3.  Robin  Boundary  Conditions. 


4  p 
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Table  4 


Maximum  allowable  CFL-number  obtained  from  direct  numerical  simulation  of  Burgers  equation. 


■Rsa 

Exact  BC 

4.50 

Penalty  BC 

6.75 

Note,  that  the  only  difference  between  the  p2irameter  values  quoted  here  and  those 
found  is  Lemma  3.2,  is  a  factor  of  1/4  on  those  terms  related  to  the  diffusion  operator. 
This  reduction  is  found  to  lead  to  optimal  time-step  restrictions. 

We  would  like  to  stress  the  importance  of  choosing  the  appropriate  value  of  the 
penalty  parameter.  It  is  our  experience,  that  this  is  best  done  by  deriving  the  theoret¬ 
ical  value  of  this  parameter  through  an  analysis  similar  to  that  done  in  Sec.  3.1.  This 
leads  to  a  parameter  which  scales  correctly  with  the  resolution  and  other  significant 
parameters.  If  the  time-step  restriction  is  dominated  by  a  viscous  time-scale,  it  is  very 
likely  that  the  theoretical  estimate  leads  to  severe  time-step  restrictions.  However, 
the  theoretical  value  may  often  be  decreased  considerably,  and  good  results  may  be 
obtained  after  only  a  few  tests.  As  we  have  seen  for  Burgers  equation,  decreasing  the 
penalty  parameter  four  times  leads  to  acceptable  Ci^L-restrictions.  We  are  not  aware 
of  any  systematic  way  of  determining  the  optimal  factor  by  which  the  theoretical  value 
should  be  decreased,  but  it  may  usually  be  determined  by  trial  and  error  through  a 
few  tests. 

To  conclude  our  study  we  have  solved  Burgers  (Eq.(3))  with  initial  condition 


(12)  £f(x,0)  =  (l-x)(l-x2)  , 

and  homogeneous  Dirichlet  boundary  conditions.  A  typical  temporal  evolution  is 
shown  in  Fig.  3.  In  Table  4  we  show  the  maximum  CFL-number  resulting  in  a  stable 
scheme.  These  results  confirm  that  the  results  from  the  linear  analysis  carries  over  to 
the  non-linear  problem. 

5.  The  Compressible  Navier-Stokes  Equations.  In  the  present  section,  we 
obtain  energy  estimates  for  the  solution  to  the  three-dimensional  compressible  Navier- 
Stokes  equations  given  in  conservation  form.  Additionally,  we  derive  open  boundary 
conditions  taking  into  account  the  full  stress-tensor,  and  prove  well-posedness  for  the 
continuous  problem.  The  derivations  follow  the  approach  introduced  in  [8,  9].  The 
main  difference  being  that  we  develop  the  theory  for  the  conservation  form  of  the 
Navier-Stokes  equations  and  that  we  include  the  off-diagonal  terms  of  the  stress-tensor 
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Fig.  3.  Temporal  evolution  of  Burgers  equation  with  initial  conditions  given  by  Eq.(li). 

in  the  full  derivations.  In  the  second  part  of  this  section,  we  continue  by  showing  how 
to  apply  the  boundary  conditions  and  prove  asymptotic  stability  of  the  semi-discrete 
scheme. 

Consider  now  the  non-dimensionalized,  compressible  Navier-Stokes  equations  given 
in  conservation  form 

aq  dF  dG  an  i  raF^  aa^  auA 

with  X  €  n  =  [-1, 1]^.  The  state  vector,  q,  and  the  inviscid  flux  vectors  are  given  as 


p 

pu 

pv 

pw 

pu 

4-p 

puv 

puw 

pv 

,  F  = 

puv 

,  G  = 

pv"^  +P 

,  H  = 

pvw 

pw 

puw 

pvw 

pw^  -f-  p 

E  . 

(^-l-p)lt 

_  {E  +  p)v 

{E  -1-  p)w 

Here  p  is  the  density,  it,  u,  w  are  the  three  Cartesian  velocity  components,  E  is  the 
total  energy  and  p  is  the  pressure.  In  the  remaining  part  of  the  paper  we  will  use 
{x,y,z)  and  (z],Z2>sb3)  interchangeably  to  denote  the  spatial  coordinates.  The  total 
energy 

E  =  p(t+^{u^  +  v^  +  w^)'^  , 
and  the  pressure  is  related  through  the  ideal  gas  law 

P  =  (7  -  l)pT  , 
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where  T  is  the  temperr»ture  field  and  7  =  Cp/c„  is  the  ratio  between  the  heat  capacities 
at  constant  pressure  (cp)  and  volume  (c„),  respectively. 

The  viscous  flux  vectors  are  given  as 


F.= 


I  = 


TxxU  +  TyxV  +  TxxW  + 


+  TyyV  +  TzyW  + 


L  rxzU  +  TyxV  +  TxzW  \ 

Considering  only  Newtonian  fluids,  the  stress  tensor  elements  are  given  as 


/  du  dv  dw\ 

fdu  dv\ 

^  dy^  dz  ) 

t  '^xy  —  Tyi 

n  9V  1 

fdu  dv  dw\ 

fdw  dv\ 

{dx  ^  dy  dz) 

>  ~  "^zy 

='‘U+fe)  ■ 

/du  dv  dw\ 

/  dw  du\ 

Vdx  ^  dy^  dz) 

t  '^xz  ~  "^zx 

where  /i  is  the  dynamic  viscosity,  A  is  the  bulk  viscosity  and  k  is  the  coefficient  of 
thermal  conductivity. 

The  equations  are  normalized  using  the  reference  values,  ttref  =  xtO)  ^ref  = 
POf  Pref  =  '^ref  =  ‘^ol^v  and  a  reference  length  L,  where  (po>'“o)  is  some  uni¬ 

form  state,  e.g.  the  ambient  free-stream  conditions  of  the  flow.  This  gives  a  Reynolds 
number  as  Re  =  pQUoLfuo  and  a  Prandtl  number  as  Pr  =  CpPo/ko. 

5.1.  Well-Posedness  and  Open  Boundary  Conditions  for  the  Contin¬ 
uous  Problem.  Consider  the  linearized,  constant  coefficient  form  of  Eq.(13).  The 
viscous  fluxes  are  split  as 


F;.  +  Fi,  +  F5 


>~dx 
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0 

0 

\  dw 

+ 

0 

0 

u— 

dz 

.  +  f*”!;  - 

G„  =  Gp  +  +  G^  - 


u— 

May 


(A  +  2/x)|| 

u— 

May 


L  M«it.+  ■'■  M“’W  +  feS 


0 

0 

mI^ 

0 

\  du 
^dx 

+ 

\  dw 
''  az 

0 

u— 

Maz 

.  ■'»S + , 

H, 


Hp  +  HJf  +  - 


tt 

M^ 

u— 

Maz 


(A  +  2/i)|^ 

M«f^+M^fl  +  (A  +  2M)«»t  +  Sf 


0 

^  0 

,,dw 

M-g^ 

0 

+ 

0 

+ 

9w 

M-gii 

\  9tJ 

.  +  M«|7  - 

Introducing  the  transformation  Jacobians 


.  aF  .  aG 
~  dq  ’  '  aq 


dFp  „  dGp 

9q, 


^33 


aHp 

aqz 


I  ^^23 


,  ^13  = 


/aF^  aH^ 

\  aq^  aqi 
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allows  for  writing  Navier-Stokes  equations  as 


at  oxi 

>=i 


,  3  3 


i=l  j~i 


d\ 

dx\  dxj 


It  is  well  known  that  Navier-Stokes  equations,  although  not  of  hyperbolic  nature,  sup¬ 
port  waves  very  similar  to  those  encountered  in  the  hyperbolic  Eider  equations.  For 
hyperbolic  systems,  Gottlieb  et  al.  [17]  have  shown  that  enforcing  the  boundary  condi¬ 
tions  through  the  characteristic  variables  of  the  system  results  a  stable  approximation. 

For  Navier-Stokes  equations,  we  linearize  around  a  uniform  state,  qo,  by  fixing  all 
the  matrices.  We  transform  into  characteristic  variables  by  diagonalizing  Ai  through 
a  similarity  transform  A  =  4S~Mi«S,  where  A  is  the  eigenvalue  matrix  and  S  and  S~^ 
are  the  matrix  of  right  and  left  eigenvectors,  respectively.  These  matrices  are  given  in 
the  Appendix.  Applying  this,  the  symmetrized,  linearized  set  of  equations  transforms 
into 


(14) 


i=l 


dxi  Re 


.=1  j=i 


dx^  dxj 


where  R  =  <5  are  the  characteristic  variables.  We  have  introduced  a  positive 
definite,  symmetrizing  diagonal  matrix,  QTQ,  given  as 


1 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

y-i 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

1 

where  cq  =  yJ'lPo! Po  is  the  uniform  state  sound  speed.  Also  we  define  the  symmetrized 
matrices 


At  =  Q^QS'UiS  ,  Btj  =  Q^QS-^Bi,S  . 


The  explicit  form  of  the  symmetric  matrices  are  given  in  the  Appendix.  The  charac¬ 
teristic  variables,  R  =  [Ri,  R-z,  R3,  R4,  Rs]^,  are  given  as 

r  pu  -  Uop  ^  (^E  +  +  v'^  +  wl)p  -  Uopu  -  Vopv  -  Wopw^  1 


R  = 


pv  -  VoP 

p-^  \p{ul  +  +  Wo)  -  “opv  -  Vopv  -  Wopw) 

pw  -  WqP 

-{pu  -  Uop)  -I-  ^  (r  -1-  \p{ul  -I-  +  w^)  -  Uopu  -  Vopv  -  Wopw') 
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We  are  now  ready  to  state  the  following 


Lemma  5.1.  Assume  there  exists  a  solution,  q,  which  is  periodic  or  held  at  a 
constant  value  at  the  y-  and  z-boundary.  If  the  boundary  conditions  in  the  x-direction 
are  given  such  that 


V(y,2:)  €  fly  X  :  -- 


dR 


J=l 


<0  , 


x=— 1 


and  the  fluid  properties  are  constrained  by 


/^>0>  ^o<0,  Ao  +  /io>0>  >  0  »  7  >  1  > 


then  Eg.  (14)  is  a  well-posed  problem  and  the  solution  is  bounded  as 

^  ^  dsi  <  0 

-r  cfxi 
=» 

Proof.  Construct  the  energy  integral  as 


Ijd 

2dt 


\mf  = 


RMJR 


Bust: 


j=l 


-ax. 


1 

dydz 

X=—\ 


tfxi 


dil 


where  0  =  fix  x  fiy  x  O^.  In  deriving  this  expression,  we  use  partial  integration 
and  assume  the  solution  to  be  periodic  or  held  at  a  constant  value  along  the  y-  and 
z-boundaries,  i.e.  contributions  from  these  boundaries  cancel.  This  is  not  a  severe 
restriction,  as  this  assumption  is  valid  for  a  large  variety  of  situations  where  open 
boundary  conditions  are  applied. 

It  is  evident  that  if  we  can  prove 


(15) 


aR^  aR 

dxi  ^'^dxj 


dn  >  0  , 


then  weU-posedness  may  be  ensured  by  properly  constructing  the  boundary  operator 
at  the  x-boundary. 
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Since  the  matrices,  are  all  symmetric,  Eq.(15)  may  be  rewritten  in  a  block- 
quadratic  form  as 


4-  / 

!Re  Jq 


WRdil  >  0  , 


where  we  have  introduced 


-  raa  dR 

^~[dx'  dy'  dz\  >  ^  - 


2BU  BU 

B\2  2^22  ^23 

B\3  B23  2^33 


We  observe  that  W*  is  a  15  X  15  symmetric  matrix,  ensuring  that  the  eigenvalue 
spectrum,  p{H*),  is  real.  Hence,  if  is  positive  semi-definite,  Eq.(15)  is  obeyed.  The 
eigenvalue  spectrum,  p(W*),  may  be  found  to  be 


Pi 

=  P2  =  P3  =  0 

P4 

0 

1 

tl 

Ps 

—  P%~  2(Ao  +  3^) 

P7 

=  Ps  =  3/io  -  \ 

//i^  -f-  2(/io  +  Ao)2 

P9 

=  Pio  =  3/io  + 

\//ig -b  2(/io  +  Ao)2 

Pli 

=  7/io  +  4Ao  — 

\J +  4(/io  +  Ao)(3/io  +  2Ao) 

Pi2 

=  7/io  +  4Ao  + 

\J 4(/io  -b  Ao)(3/io  -b  2Ao) 

(  2cg  .  2(7  -  l)*:o 

Pl3 

=  Pl4  =  Pih  =  \ 

Pr  • 

Here  subscript  ’0’  signifies  the  parameter  values  in  the  uniform  state  around  which  we 
have  linearized.  For  most  real  fluids  under  non-extreme  conditions,  it  is  true  that  /xo 
is  positive,  Aq  is  negative  and  the  following  relationship  is  obeyed  [18] 


(16)  >  Ao  -b  2fM}  >  fM)  ■ 

A  simple  investigation  of  the  eigenvalues  reveals  that  W  is  positive  semi-definite  under 
these  conditions.  Thus,  Eq.(15)  is  true  provided 

Po  ^  0  ,  Ao<0,  Ao-b/io^Oi  —  6  ,  7^1- 

Pr 

These  conditions  are  only  natural  as  discussed  in  [19].  In  fact,  if  they  are  not  obeyed, 
Navier-Stokes  equations  violates  the  second  law  of  thermodynamics. 
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We  now  obtain  that  well-posedness  is  ensured  under  the  additional  condition 

<0  , 


v(y,  z)  e  rij,  X  -  - 


and  the  solution  is  bounded  as 


„  2  -r  5R 


J  =  1 


X=:  — 1 


1-1 

2di 


~]dn<0  , 


where  QR  =  ’q.  □ 


As  stated  in  Lemma  5.1,  appropriate  boundary  conditions  at  the  x-boundary  have 
to  obey 

T  1 

<  0  . 


„  2  ^  dR 


We  now  define 


i=i 


j=l 


-1 


where 


Q  ^0(7  ~  ^)  ^Ci  ■^0  +  2/io  8^2  .  Aq  +  /xq  / dR2  dRj  \ 
’  2poPr  dx  2po  dx  po  \  dy  dz  J 


Po  dR2  ,  Ao  +  A*o  9(2 
tr2  = - ^  + 


po  dx 


Po  dy 


(17) 


G3  =  - 

G4  = 


fco(7  -  1)  dCi 

2poCoPr 

/io  dR^  ^  Ao  +  /io  ^2 


Po  dx  '  Po  dz 

Q  _feo(7~l)^Ci  Aq  +  2/xo  ^^2  yl•o  +  p^o^dR2 

*  2poPr  dx  2po  dx  po  \  dy  dz  )  ' 

where  we,  for  simplicity,  have  introduced 


Cl  =  fZj  +  iZs  - 


2co 

7^ 


R3  i  C2  =  R\  —  R5 


This  allows  lor  rewriting  the  constraint  on  the  boundary  contribution  as 


R^AR-  -^R^G 
Re 


Q  <  0  , 


-1 


where  A  is  the  diagonal  eigenvalue  matrix  obtained  from  the  similarity  transform.  We 
now  reformulate  this  as 

('«)  EV((lW-£l^G,)'-(£G0')  <0, 
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where  are  the  wave  speeds  by  which  the  characteristic  variables  are  advected,  as 
given  by  the  diagonal  elements  of  A,  cind  we  have  introduced  e  =  Re“^ ,  This  formula¬ 
tion  makes  it  straightforward  to  devise  inflow-outflow  boundary  conditions,  which  axe 
maximal  dissipative  and  ensure  weU-posedness  of  the  complete  problem. 

We  note  in  particular  that  this  formulation  takes  into  account  the  off-diagonal 
terms  of  the  stress  tensor,  which  is  neglected  in  most  previous  work  [8,  9,  10].  These 
terms  may  be  of  importance  if  the  artiflcial  boundary  is  introduced  into  a  strongly 
vortical  region  of  the  flow,  e.g.  a  wake  flow  behind  a  blunt  body. 

Inflow  Boundary  Conditions.  At  x  =  —1,  Eq.(18)  becomes 

§  E  ^r'  ( <  0  • 

Subsonic  Inflow  :  A^  >  0  ,  Aa  >  0  ,  A3  >  0  ,  A4  >  0  ,  As  <  0 

(19)  AiiZi  —  eG\  =  0 

^2^2  ~  —  0 

A3R3  —  eGs  =  0 

X4R4  —  £(?4  =  0 

eGs  =  0 

Supersonic  Inflow  :  Ai  >  0  ,  A2  >  0  ,  A3  >  0  ,  A4  >  0  ,  As  >  0 

(20)  AifZj  —  eGi  =  0 

^2-^2  ~  £^2  =  0 
^3^3  ~  £^3  =  0 
A4R4  —  £^4  =  0 
AsiZs  ~  £Gs  =  0 

Outflow  Boundary  Conditions.  At  x  =  1,  Eq.(18)  becomes 

I E  ^  0  ■ 

Subsonic  Outflow  :  Aj  >  0  ,  A2  >  0  ,  A3  >  0  ,  A4  >  0  ,  As  <  0 

(21)  eG2  =  0 

eGs  =  0 

£^4  =  0 

(AsjRs  +  £^5  =  0 
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Supersonic  Outflow  :  Ai  >  0  ,  A2  >  0  ,  A3  >  0  ,  A4  >  0  ,  A5  >  0 

eGa  =  0 

cGs  =  0 
or 

CG4  =  0 
eGs  =  0 

We  note  that  for  both  types  of  outflow  boundary  conditions,  it  is  only  necessary  to 
specify  four  conditions,  since  eGa  =  0  eGi  =  —eGs.  Due  to  the  special  structure 
of  G  we  also  observe  that  adding  an  extra  condition  on  eGi  does  not  place  extra 
conditions  on  the  solution,  since  such  a  condition  is  redundant.  This  observation  will 
be  used  later. 

It  was  shown  by  Strikwerda  [20]  that  the  proper  number  of  boundary  conditions 
for  an  incomplete,  parabolic  system,  like  the  compressible  Navier-Stokes  equations,  is 
5  in  the  inflow  region  and  4  in  the  outflow  region.  Our  result  clearly  conforms  with 
that. 

We  cdso  note  that  in  the  limit  of  infinite  Reynolds  number,  these  boundary  con¬ 
ditions  converge  uniformly  toward  the  well  known  characteristic  boundary  conditions 
for  the  compressible  Euler  equations  [21].  This  property  is  important  in  order  to  avoid 
weak  boundary  layers  of  the  order  exp(—x/e)  (see  [8]). 

5.2.  The  Semi-Discrete  Scheme.  Following  the  line  of  thought  that  led  to  the 
asymptotically  stable  scheme  for  Burgers  equation,  we  propose  a  Legendre  collocation 
scheme  for  enforcing  open  boundary  conditions  to  the  compressible  Navier-Stokes 
equations 

^  ^  dx'*'  dy  dz  ~  Re\dx  dy  dz  J 

-r,Q-(x)S  {n-  (R  -  S-^g4(t))  - 

-r2Q^ix)s  (r  -  s-^g2(t))  +  • 

Here  Q~(x)  and  are  given  by  Eq.(9)  and  S  is  the  right  eigenvector  matrix  as 

given  in  the  Appendix.  The  boundary  conditions  for  the  state  vector  are  given  through 
the  two  vectors,  gi(t)  and  g2(t),  which  we  for  convenience  assume  to  be  uniform.  The 
four  matrices,  TZ~,  Q~  and  are  chosen  such  as  to  construct  the  appropriate 
boundary  operator  as  derived  in  the  previous  section.  Hence,  we  have  for  the  inflow 
region 


(22) 


sGi  =  0 
tGi  —  0 
tGz  =  0 
CG4  =  0 
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•^1 

0 

0 

0 

0 

■  1 

0 

0 

0 

0  ■ 

0 

A2 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

■^3 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

A4 

0 

0 

0 

0 

1 

0 

,  0 

0 

0 

0 

aAs 

0 

0 

0 

0 

1 

where  a  =  0  for  subsonic  conditions  and  a  =  1  for  supersonic  conditions.  Likewise  we 
define 


71+  = 


o 

0 

0 

0 

o 

o 

0 

0 

0 

f . 

o 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

.  0-^  = 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

o 

0 

0 

0 

.  0 

0 

0 

0 

1 . 
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where  =  1  for  subsonic  conditions  and  /?  =  0  for  supersonic  outflow  conditions. 

We  have  to  choose  rj  and  T2  such  that  the  semi-discrete  scheme  is  asymptotically 
stable.  The  proper  choice  is  stated  in  the  following  Lemma. 


Lemma  5.2.  Assume  there  exists  a  solution,  q,  which  is  periodic  or  held  at  a 
constant  value  at  the  y-  and  z-boundary,  and  that  the  fluid  properties  of  the  uniform 
state,  qo,  are  constrained  by 

Po>0,  Ao<0,  Ao  +  /io>0.  >  0  .  7  >  1  > 

and  related  as 


>  •^o  +  2/io  >  fM)  • 

The  linearized,  constant  coefficient  version  of  the  scheme  given  by  Eq.(23)  is  asymp¬ 
totically  stable  at  the  inflow  if 

(l  -I-  A  -I-  Vl  +  «)  >  Ti  >  —  (l  -I-  K  -  Vl  +  k)  • 

Here 

e  ko 

K  = - . 

2u  Prpo'Uo 

These  results  are  independent  on  whether  the  inflow  is  subsonic  or  supersonic. 
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For  supersonic  outflow 


For  subsonic  outflow 


The  solution  to  Eq.(23)  is  hounded  in  the  form 


1.  J  J  x~xit 

Proof  Write  Eq.(23)  is  its  symmetrized,  linearized,  constant  coefficient  version 


2  2 at 


-r,(3'(x)  (Q^Q7^-R- 
-T2Q-^ix)(QFQn+K+~Q^QQ+G^  , 


where  we,  without  loss  of  generality,  have  assumed  homogeneous  boundary  conditions. 
We  construct  the  energy  integral,  apply  the  Gauss- Lobatto  quadrature  rule  and  partial 
integration  to  obtain 


(24|i||QR||5,  =  I  -1 

-Tiwjf  K^Q^QIZ-R-eK^g-'p^Bl—  dydz 

-Tiojj  j  R^Q^Qn+K  +  eR^g+Y.^h^  ^ 

**  .  ^  J  X=1 


where  we  have  used  the  assumption  about  periodicity  or  constant  value  at  the  y-  and 
z-boundary.  Additionally,  we  have  introduced  e  =  Re~'  and  oj,  which  is  the  Legendre 
weight  at  the  endpoints  and  applied  the  definition 

j=x  J 
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Using  the  Gauss-Lobatto  quadrature  rule  allows  for  writing 


dO.  = 


/  / 

JUyJu. 


^  ^  dK^  ,  dR 

dxi  '^dxj 
•=i  j=»  J 


X=Xk 


+  w 


^  ^  ^5: 

dxi  ’^ax, 

»=1  ]—X  ■> 

^  ^  K* 

4^^  ax.  ’^ax, 

t=l  j=t  *  J 


ars— 1 


xsl  > 


dydz  >0  . 


Here  x^  signifies  the  Legendre- Gauss-Lobatto  collocation  points.  The  inequality  fol¬ 
lows  from  the  analysis  done  in  the  proof  for  Lemma  5.1,  and  is  ensured  provided  the 
fiuid  properties  are  constrained  by 

/xo>0,  Ao<0,  Ao  +  /xo>0,  ^>0,  7>1- 

Pr 

It  was  shown  by  Abarbanel  and  Gottlieb  [18]  that  if  a  scheme  is  stable  without  the 
contributions  from  the  off-diagonal  stress- tensor  terms,  then  it  will  rer  ain  so  even  if 
the  these  terms  are  included.  This  is  a  consequence  of  the  general  relation 


>  Ao  4-  2/io  >  IM)  , 

which  roughly  gives  the  relation  between  the  eigenvalues  of  the  normal  stress-tensor 
elements  and  the  off-diagonal  elements.  Thus,  it  is  sufficient  to  prove  stability  in  the 
absence  of  the  off-diagonal  contributions. 

The  penalty  parameters,  ri  and  T2,  has  to  be  chosen  such  that  the  boundary  term 
of  the  energy  integral  not  destroys  the  stability  of  the  Cauchy-problem.  We  treat  the 
two  boundary  contributions  separately. 


Inflow  Condition.  The  contribution  of  the  boundary  term  at  the  inflow  (x  =  -1)  fol¬ 
lows  from  combining  Eq.(24)  and  Eq.(25)  and  neglecting  the  off-diagonal  contributions 
to  obtain 


R’’  R  -  £R^  (I  -  nwC')  s;,  ^  <  o  , 


where  I  is  the  identity  matrix. 
First  we  note  that 

dR^ 


-tU)- 


B: 


9X2  ^^9x2 


9R  dR 
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since  B22  ^^33  positive  semi-deiinite  with  an  eigenvalue  spectrum  given  as 


=  PliBis)  =  0  p2iB^22)  =  P2m3)  =  ^ 

P3{Bh)  =  P3{B^3s)  =  2^  P4iBh)  =  P4{Bi:,)  =  2^  . 

^5(^22)  =  +  1)  PsiBis)  =  +  1) 

Since  all  matrices  are  symmetric,  the  remaining  part  of  the  constraint  may  be  expressed 
in  block-quadratic  form  as 


R^n-K  <  0  , 


where 


R  = 

«  9R 

-s- 

T 

Ai  —  2t\uQ^ QTZ  -e(l  -  Txu)B\.^ 

’  dx 

2 

— e(l  -  Tju)Bfi  -2euBii 

where  we  have  used  Q~  =  I.  H~  is  a  10  x  10  symmetric  block-matrix.  Similar  to 
the  approach  applied  in  Sec.  3.2,  we  have  transformed  the  problem  of  stability  into 
proving  that  'H~,  for  a  suitable  value  of  Tj,  is  negative  semi-definite.  The  eigenvalue- 
spectrum,  p(H~),  can  be  found  by  doing  a  LU-decomposition.  Since  H~  is  symmetric, 
the  eigenvalues  appear  as  pi{'H~)  =  Ua- 

We  will  not  give  the  general  form  of  the  eigenvalues  here,  since  they  are  rather 
complicated.  However,  straightforward  but  very  lengthy  algebra  shows  that  all  eigen¬ 
values  are  negative  if  t\  is  chosen  such  that 


(1  -1-  K  +  \/l  +  k)  >  n  >  —  (1  +  K  -  -v/l  +  «)  , 

UK  '  ■'  UK  ^  ' 


where 


-  _L  ^0 

2u  Fipouo 

This  result  is  independent  of  whether  the  inflow  is  subsonic  or  supersonic. 


Outflow  Condition.  Neglecting  the  contribution  from  the  off-diagonal  terms  yields  a 
criteria  for  stability  at  the  outflow  (x  =  1) 

-a’'  +  T2U.C’'QK+)  R  +  eR'’  (I  -  T,u-ff+)  E  s  0  • 

Similar  to  the  approach  followed  in  the  previous  part  of  the  proof,  we  see  that  the 
contributions  from  B22  and  B^^  are  always  negative  and  independently  ensure  stability. 

We  now  rewrite  the  remaining  part  of  the  condition  at  the  outflow  in  block- 
quadratic  form; 

R^?f+R  <  0  , 
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where 


n+  =  -  ~ 

2  £(1  -  Tiu)B\i  -2eufBii 

To  form  we  have  assumed  =  T.  The.  additional  boundary  condition  introduced 
by  this  replacement  is  redundant  as  discussed  in  Sec.  5.1.2.,  and,  hence,  no  extra 
restrictions  are  put  on  the  system  by  this  approach.  The  eigenvalue  spectrum, 
may  again  be  found  through  a  LU-decomposition.  We  state  here  only  the  bounds  on 
T2  that  ensure  negative  semi-definiteness  of  for  supersonic  outflow 


For  subsonic  outflow  the  bounds  become 


Combining  £q.(24)  and  Eq.(25),  we  obtain  a  bound  for  the  growth  of  the  solution 


We  wish  to  emphasize  that  the  bounds  on  Ti  and  given  in  Lemma  5.2  remains 
valid  in  the  limit  when  the  Reynolds  number  approaches  infinity.  This  is  easily  realized 
by  expanding  the  bounds  for  c  <  1  to  obtain 


^  1  1 

OO  >  Ti  >  —  +S—K  , 

2w  8w  ’ 


in  the  inflow  region  and 


OO  >  T2  >  -OO 


OO  >  T2  > 


2b; 
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for  supersonic  and  subsonic  outflow,  respectively.  The  linewzed,  constant  coeflicient 
version  of  the  Euler  equations  may  be  transformed  into  5  independent  hyperbolic  equa¬ 
tions  for  which  we  should  expect  the  bounds  on  the  penalty  parameters  to  be  given 
by  the  restdts  in  Sec.  3.1.1.  We  observe  that  the  bounds  given  above  converge  uni¬ 
formly  to  the  expected  values  in  the  limit  of  vanishing  viscosity  and,  thus,  the  scheme 
remains  stable.  The  observation  that  no  bounds  are  necessary  on  T2  for  supersonic 
outflow  simply  reflects  the  fact  that  no  boimdarj  conditions  are  required  for  the  Euler 
equations  at  such  a  boundary. 
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5.3.  Numerical  Tests.  The  proof  given  in  the  previous  section  is  only  strictly 
valid  for  the  linearized,  constant  coefficient  version  of  Navier-Stokes  equations.  To 
vahdate  the  results  and  show  that  it  carries  over  to  the  full  non-linear  Navier-Stokes 
equations,  we  have  implemented  the  scheme  in  an  existing  spectral  code  (see  [22]  for 
details),  originedly  developed  for  studying  two-dimensional  compressible  flow  around 
an  infinitely  long  circular  cylinder. 

As  spatial  approximation  scheme  was  used  a  standard  Fourier- Chebyshev  collo¬ 
cation  scheme  in  polar  coordinates,  (r,  6),  with  a  3rd-order  Runge-Kutta  method  for 
time-stepping. 

The  new  scheme  is  simple  to  implement  in  existing  codes,  as  we  only  need  to  apply 
a  correction  of  the  flux  of  the  state  vector  at  the  boundary.  Following  the  scheme, 
given  by  Eq.(23),  we  need  to  derive  the  two  vectors  R  and  G.  The  characteristic 
variables  are  given  as 


Ri  =  (m,  -  pur)  -I-  —  , 

Co 

^2  =  P  -  ^  , 

Rz  =  me  —  pug  , 

R4  =  -{mr-pUr)  +  —  . 

Co 

where  cq  is  the  uniform  state  sound  speed. 

We  have  for  convenience  introduced 

Ur  =  Uoki  +  Vok2  ,  Ug  =  140*2  “  ^0*1  , 

which  are  the  radial  and  azimuthal  velocity  components,  respectively,  of  the  uniform 
state  and 

TTtr  “  ITlfiki  771^*2  »  TTlg  —  771,^ *2  —  771^*1  , 

are  the  radial  and  azimuthal  components  of  the  momentum  of  the  flow  field.  Here  k  = 
(*i)  *2)  signifies  an  outward  pointing  normal- vector  at  the  boundary.  The  linearized 
pressure,  p,  is  given  as 

P  =  (7  -  1)  -S  +  ^P(«o  +  »o)  -  ■«*oT7i„  -  Vom„  . 

The  eigenvalues  corresponding  to  the  characteristic  functions  and  determining  the 
direction  and  propagation  velocity  of  the  characteristic  waves,  are 

Ai  =  Ur  -J-  Co  ,  A2  =  A3  =  lir  ,  A4  =  Ur  —  Co  . 
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Following  the  approach  outlined  in  the  previous  section,  we  have  likewise  derived  the 
viscous  correction  vector,  G,  at  the  outer  boundary  as 


Q  _  (7 -1)^0  1  dCi  2  ftp  8(2  IfiodRs 

^  Pr  2po  dr  3  po  dr  3  po  dd  ’ 

(7  -  1)^0  1  dC\ 

^  Pr  2cqPq  dr  ’ 

f-i  _  tM)  dRs  1  po  d(2 

Po  dr  6  Po  dd  ’ 

^(7-l)ko  1  aCi  2fiodC2  jlpodRa 

Pr  2po  dr  3  po  dr  3  po  dd  ’ 

where  again  we  have  defined 


C\  ~  R\  ■\-  Ra - ^^R2  »  C2  —  R\—Ra  • 

7-1 

Also,  we  have  J  =  1/r,  which  is  the  transformation  Jacobian  from  Cartesian  to  polar 
coordinates.  We  note  that  no  extra  calculation  of  derivatives  is  needed  in  order  to 
form  the  two  vectors,  since  the  radial  and  azimuthal  derivatives  at  the  boundary  are 
calculated  during  evaluation  of  the  interior  dynamics  when  employing  a  global  scheme. 
Thus,  the  only  additional  requirement  is  to  store  values  of  the  derivatives  of  the  state 
vector  at  the  boundary,  i.e.  the  computational  requirement  for  enforcing  this  new 
method  is  negligible. 

The  boundary  conditions  are  enforced  at  each  intermediate  time  step  of  the  Runge- 
Kutta  method.  Simulations  were  done  with  a  Reynolds  number  of  100,  a  Mach  number 
of  0.4,  the  diameter  {D)  of  the  cylinder  being  6.10  cm  and  the  reference  temperature 
was  300° AT.  These  parameters  ensure  that  the  flow  field  remains  subsonic.  The  reso¬ 
lution  was  96  Fourier-modes,  72  Chebyshev  modes  and  the  radius  (£)  of  the  compu¬ 
tational  domain  was  20  cylinder  diameters. 

As  penalty  parameters  we  used 


AT*  2,, 


2 


vnn;) ,  t3  =  — -  , 


where  N  is  the  number  of  Chebyshev  modes,  2/Ir  is  a  factor  occurring  from  the  radial 
mapping  of  L  into  [-1,1]  and 

_  eN^  kp 

2  PipoUo  ’ 

This  choice  appears  naturally  from  the  results  stated  in  Lemma  5.2,  and  the 
experience  gained  in  Sec.  4,  indicating  that  for  dissipative  terms  we  should  reduce  by 
a  factor  of  4  in  order  to  obtain  the  optimal  value  of  r.  With  this  choice  of  penalty 


34 


parameters  we  were  able  to  perform  the  simulations  without  any  reduction  in  time- 
step  as  compared  to  the  exact  method  of  enforcing  the  boundary  conditions.  It  shoxild 
be  mentioned,  that  in  the  original  code  only  characteristic  boundary  conditions  for 
the  Euler  equations  were  enforced.  Comparing  with  results  discussed  discussed  in  Sec. 
4,  we  observe  that  for  3rd-order  Runge-Kutta  we  should  expect  the  two  methods  to 
impose  almost  equivalent  time-step  restrictions.  This  is  confirmed  by  the  simulations 
<ind  shows  that  the  results  from  the  simple  linear  analysis  carries  over  to  the  full 
non-linear  Navier-Stokes  equations  in  this  case. 

In  Fig.  4  we  show  contour-plots  of  the  normalized  density  and  the  pressure  at 
T=143.5,  corresponding  to  approximately  23  shedding  cycles.  The  von  Karman  vortex 
street  is  clearly  demonstrated,  and  we  observe  that  the  boundary  conditions  at  the 
outflow  boundary  affect  the  flow  only  slightly.  The  Strouhal  number  for  the  shredding 
frequency  is  found  to  be  St  =  0.163,  which  is  in  full  accordance  with  experimental 
findings  [23]  and  we  observe  no  spurious  frequencies  or  reflections  from  the  artificial 
boundary  back  into  the  flow  field  (see  [22]  for  a  further  discussion  of  this). 


Pig.  4.  Contour  plots  of  the  normalized  density,  p/po,  and  the  normalized  pressure,  p/po,  at  the 
non-dimensional  time  T=14S.5  for  a  flavi  atRe  =  100,  M  =  0.4,  Z?  =  8.10  cm  and  To  =  SOO'AT. 

6.  Concluding  Remarks.  The  purpose  of  the  present  paper  has  been  two-fold. 
The  first  goal  has  been  to  develop  boundary  conditions  for  wave-dominated  problems, 
leading  to  well-posed  total  problems.  It  was  argued,  that  for  smooth  solutions  and  the 
kind  of  operators  we  have  considered  here,  it  is  sufficient  to  consider  the  problem  of 
well-posedness  for  the  linearized,  constant  coefficient  version  of  the  non-linear  initiaJ- 
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boundary  value  problem.  Using  this  allowed  for  deriving  proper  boundary  conditions  to  Burgers 
equation  and  the  three-dimensional,  compressible  Navier-Stokes  equations,  and  these  boundary 
conditions  were  shown  to  ensure  well-posedness  of  the  total  problem.  It  should  be  stressed  that 
the  boundary  conditions  derived  for  the  Navier-Stokes  equations  takes  into  account  all  elements  of 
the  stress-tensor,  and  only  very  light  assumptions  were  made  to  derive  these.  Additionally,  they 
remain  valid  even  in  the  limit  of  vanishing  viscosity. 

Having  derived  appropriate  boundary  conditions  naturally  leads  to  the  question  of  how  to 
enforce  these  in  a  discrete  approximation  of  the  problem.  This  hats  been  the  second,  and  main, 
contribution  of  the  paper.  Recent  results  [7]  on  the  connection  between  stability  of  discrete  and 
semi-discrete  approximations,  suggest  that  it  is  sufficient  to  consider  asymptotic  stability  for  the 
semi-discrete  approximation.  We  have  only  considered  Legendre  collocation  methods  here.  This 
choice  is  merely  dictated  by  a  wish  to  obtain  analytical  results  and  we  have  indicated,  by  numerical 
tests,  that  all  results  carry  over  to  Chebyshev  collocation  operators.  The  stability  proofs  for  the 
semi-discrete  approximations  to  the  linearized,  constant  coefficient  versions  of  Burgers  equation  and 
the  compressible  Navier-Stokes  equations  are  all  completed  by  using  the  classical  energy  method. 
We  emphasized  that  the  proposed  schemes  remain  stable  even  in  the  limit  where  the  problems 
become  purely  hyperbolic. 

The  proposed  penalty  method  changes  the  eigenvalue  spectra  of  the  discrete  approximations 
of  the  operators  considerably.  In  order  to  understand  this,  we  performed  a  detailed  investigation 
of  the  effect  on  the  eigenvalue  spectra  of  linear  operators.  It  has  been  shown  that  the  value  of 
the  penalty  parameter,  which  is  obtained  form  the  theoretical  analysis,  often  implies  that  the 
maximum  allowable  time-step  compares  unfavorable  with  that  allowed  through  more  traditional 
methods.  However,  we  discussed  in  detail  how  to  remedy  this  and  showed  that  choosing  the 
penalty  parameter  properly  may  allow  for  increasing  the  maximum  time-step  with  as  much  as  50%. 
Although  we  are  not  aware  of  a  systematic  way  of  determining  the  optimal  value  of  the  penalty 
parameter,  we  do  not  see  that  as  any  significant  disadvantage.  Our  experience  tells  that  once  the 
theoretical  values  of  the  penalty  parameters  are  obtained,  only  a  few  tests  are  needed  to  obtain  the 
optimal  value.  Additionally,  this  only  has  to  be  done  once,  and  since  only  a  few  hundred  time-steps 
are  required  to  test  whether  the  scheme  is  stable  or  not,  we  consider  this  an  insignificant  problem. 

Most  of  the  theoretical  results,  obtained  for  linearized,  constant  coefficient  versions  of  the  equa¬ 
tions,  are  confirmed  by  numerical  simulations  of  the  full  non-linear  equations.  It  is  stressed  that 
the  proposed  penalty  method  is  very  easy  to  implement  in  existing  codes,  which  is  an  attractive 
feature. 
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Although  all  results  and  nuiuerical  simulations  in  this  paper  are  obtained  using 
spectral  collocation  methods,  the  main  conclusions  carry  over  to  finite  difference/finite 
element  methods.  The  derivation  of  the  proper  boundary  operators,  be  that  for  Burg¬ 
ers  equation  or  for  the  compressible,  Navier-Stokes  equations,  is  obviously  unaffected 
by  the  choice  of  spatial  approximation  method.  The  proposed  penalty  method  for 
enforcing  the  boundary  conditions  may  be  applied  in  exactly  the  same  manner  as 
discussed  here,  when  using  alternative  spatial  discretization  methods.  The  only  dif¬ 
ference  is  the  value  of  the  penalty  parameter,  which  will  depend  strongly  on  the  order 
of  the  method.  Thus,  applying  an  other  method  requires  one  to  derive  this  penalty 
parameter.  This  may  be  done  by  an  approach  equivalent  to  the  one  utilized  here. 
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Dettori  from  Brown  University  for  many  useful  discussions. 


Appendix:  Symmetric  Matrices  for  the  Navier-Stokes  Equations..  Con¬ 
sider  the  linearized,  constant  coefficient  compressible  Navier-Stokes  equations  in  con¬ 
servation  form  given  as 


—  +  ^A  — 
dt  ‘di. 


«=i 


.=1  j=i 


dxf  dtSj 


The  matrix,  Ai,  diagonalizes  under  the  similarity  transform,  A  =  where  the 

right  eigenvector  matrix,  «S,  and  the  left  eigenvector  matrix,  5“'  ,  are  given  as 
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Here 
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Introducing  this  transformation  into  the  Navier-Stokes  equations  yields 


t=l 


dxi  Re 


»=1  j=i 


'  dxi  dxj 
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where  R  are  the  characteristic  variables  and  is  a  positive  definite,  symmetrizing 
diagonal  matrix. 

The  symmetrized  matrices 

At  =  Q^QS-^AiS  ,  =  Q^QS-^BijS 

are  given  as 
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We  have  for  convenience  introduced 


$  = 


7  -  1 
7  Pr 
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0  0  0  0 

-10  0  0 
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